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SUMMARY 


(1)  TECHNICAL  PROBLEM 

We  are  Investigating  mathematical  models  of  detonation  In  dense  fluids  and 
In  solids.  We  wish  to  obtain  a  fundamental  understanding  of  the  Interaction  of 
pressure  waves,  or  shockwaves,  with  high  explosives. 

(2)  GENERAL  METHODOLOGY 

We  are  solving  the  steady-state  dynamical  equations  of  continuum  mechan¬ 
ics,  Including  viscosity,  conductivity,  and  chemical  reactions,  to  simulate 
detonation  in  fluids  capable  of  undergoing  fast  exothermic  reactions.  We  are 
solving  the  coupled  ordinary  differential  equations  of  motion  of  molecular 
dynamics,  for  solids,  to  simulate  the  initiation  process  preceding  detonation. 
We  are  examining  simple  models  designed  to  describe  the  results  of  these  macro¬ 
scopic  and  microscopic  approaches. 

(3)  TECHNICAL  RESULTS 

We  have  obtained  detonation  wave  profiles,  using  a  realistic  dense-fluid 
equation  of  state.  We  have  compared  these  profiles  to  the  predictions  of  the 
simplified  Zeldovlch-Von  Neumann-Doerlng  model. 

We  have  developed  methods  for  simulating  the  rapid  uniaxial  compression  of 
solids  suited  to  molecular  dynamics  simulation. 

We  have  formulated  the  interatomic  distribution  function  in  solids,  for 
comparison  with  atomistic  simulations  and  use  in  kinetic  reaction  models. 

(4)  IMPLICATION  FOR  FURTHER  RESEARCH 

The  relatively  long  vibrational  relaxation  times  occurring  in  shocked 
diatomic  molecules  with  realistic  values  of  viscosity  and  thermal  conductivity 
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suggest  that  the  nonequilibrium  distribution  of  energy  Is  more  Important  than 
the  Interaction  of  ordinary  transport  properties  with  chemical  reactions.  We 
therefore  Intend  to  concentrate.  In  the  next  two  years,  on  reaction  initiation 
and  Intramolecular  energy  transfers  In  reacting  molecules  in  the  solid  phase. 


I.  INTRODUCTION 


Shock  and  detonation  waves  Involve  atomic-scale  processes  which  are 
Imperfectly  understood  and  hard  to  measure.  For  this  reason  we  have  embarked 
on  a  program,  using  modern  computational  methods,  of  studying  the  processes 
through  which  mechanical  shockwave  energy  stimulates  the  release  of  chemical 
energy  within  a  detonation  wave. 

Our  work  entails  the  correlation  of  two  points  of  view:  continuum 
mechanics  and  atomic  mechanics.  We  seek  to  understand  the  limitations  of  the 
continuum  approach  In  dealing  with  at.omic-scale  problems  and  to  develop  methods 
for  describing  these  small-scale  processes  in  a  realistic  way. 

Considerable  information  on  the  equilibrium  and  transport  properties  of 
simple  systems  Is  now  ava11able(l ,2) .  Computer  simulations  of  dense  fluids, 
coupled  with  thermodynamic  perturbation  theories(3)  based  on  simple  hard  and 
soft  sphere  model's,  make  it  possible  to  describe  the  pressure-volume- 
temperature  equation  of  state  accurately,  for  sufficiently  simple  materials. 

For  some  "realistic"  potentials,  such  as  the  Lennard-Jones  6-12  interatomic 
potential,  comprehensive  computer  simulation  studies  have  established  the 
density  and  temperature  dependence  of  the  viscosities  and  the  thermal  conducti¬ 
vity  as  wel 1(2).  Much  less  is  known  about  the  contribution  of  Intramolecular 
degrees  of  freedom  to  the  transport  properties  and  equation  of  state.  Studies 
on  molecules  like  methane  and  benzene  indicate  that  the  computer  simulation 
approach  Is  now  feasible,  even  for  molecules  containing  several  atoms(4). 

A  numerical  study,  comparing  the  continuum  and  atomistic  versions  of  a 
very  strong  dense-fluid  shockwave  made  use  of  both  the  equilibrium  and  non¬ 
equilibrium  constitutive  properties  of  the  6-12  potential,  and  indicated  that 
useful  conclusions  could  be  drawn  by  comparing  the  two  approaches (5).  A 
natural  extension  of  that  work  Is  to  consider  the  possibility  of  coupling 
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chemical  reactions  with  the  dissipative  processes  of  viscosity  and  thermal 
conductivity  which  appear  In  shockwaves.  That  coupled  combination  Is  a  detona¬ 
tion  wave. 

The  simplest  prototype  condensed-phase  high  explosive  Is  probably  nitric 
oxide,  NO.  This  molecule  forms  molecular  nitrogen  and  oxygen  In  an  exothermic 
reaction  generating  approximately  100  kllobars  pressure(6).  We  set  out  with 
the  go-1  of  simulating  the  detonation  of  liquid  NO  using  both  the  continuum  and 
atomistic  approaches.  The  continuum  calculations  are  complete.  Using  an 
equation  of  state  which  describes  the  Lennard-Jones  potential  throughout  the 
fluid  regions  of  the  phase  diagram,  and  transport  coefficients  from  the  Enskog 
theory,  we  have  generated  detonation  profiles  as  functions  of  *he  energy 
release  and  activation  energy  of  the  dense-fluid  detonation.  These  calcula¬ 
tions,  including  viscosities  and  thermal  conductivity,  are  compared  to  the 
Zeldovich-Von  Neumann-Doering  model  in  Section  II. 

Shortly  after  this  work  was  begun,  a  research  program  was  initiated  at  Los 
Alamos,  designed  to  measure  and  to  model  the  detonation  of  nitric  oxide.  This 
program  Includes  experimental,  quantum  mechanical,  and  dynamical  groups  working 
toward  a  common  goal.  Related  Los  Alamos  research,  carried  out  by  Brad  Hoi lan, 
Galen  Straub,  and  their  co-workers,  has  established  that  the  equilibration  time 
between  vibrational  and  translational  temperatures  In  dense  diatomic  fluids 
(nitrogen)  is  relatively  long(7),  just  as  in  the  gas  phase.  These  long  equili¬ 
bration  times,  of  order  100  vibration  times,  suggest  that  the  complete  combus¬ 
tion  process  is  impractically  long  for  an  atomistic  simulation.  For  this 
reason  we  have  concentrated  our  more  recent  vrerk  on  solid-phase  Initiation  and 
mode  equilibrium.  We  are  developing  solid-phase  calculations  In  which 
chemically  active  molecules  can  be  studied  under  uniaxial  and  homogeneous 


adiabatic  compression.  That  work  Is  described  In  Sections  III  and  IV.  We 
expect  to  devote  the  next  two  years'  effort  to  enhancing  our  understanding  of 
solid-phase  detonation  Initiation  from  the  atomistic  mechanical  point  of  view. 
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II.  DENSE— FLUID  DETONATION  WAVE  STRUCTURE 

The  classical  steady  one-dimensional  ZND  (for  Zeldovlch,  Von  Neumann,  and 
Doering)  model  of  detonation-wave  structure  can  be  summarized  with  the  help  of 
Figure  1.  It  Is  assumed  that  the  reactants  can  assume  the  equilibrium  unreact¬ 
ed  thermodynamic  states  shown  on  the  lower  Hugoniot  curve  as  the  result  of 
shock  compression.  Each  of  these  states  satisfies  the  condition  that  its  mass, 
momentum,  and  energy  fluxes  agree  with  those  characterizing  the  initial  state 
(state  0).  This  unreacted  Hugoniot  is  hypothetical  because  chemical  reactions 
will  gradually  transform  the  reactants,  at  least  for  strong  enough  shockwaves, 
to  hot  products. 

The  fully  reacted  products  can  also  be  described  by  a  Hugoniot  curve 
satisfying  the  conservation  laws.  This  Is  the  upper  Hugoniot.  Its  position, 
relative  to  the  lower  Hugoniot,  depends  upon  the  heat  of  reaction  Q.  The 
Zeldovlch-Von  Neumann-Doerlng  detonation  wave  model  pictures  the  detonation 
process  as  two  consecutive  steps:  first,  the  reactants  are  shock-compressed  to 
the  (Von  Neumann  spike)  state  S.  Next,  these  hot  reactants  slowly  undergo 
chemical  reactions,  proceeding  to  the  final  (Chapman-Jouguet)  state  CJ.  During 
the  reaction  phase  the  state  point  follows  the  Rayleigh  line,  satisfying  the 
constraints  of  constant  mass  and  momentum  flux  necessary  In  a  steady  wave. 
Neither  the  transport  coefficients  nor  the  reaction  rate  enter  into  Figure  1. 
These  state  functions  determine  the  space  and  time  scales  associated  with  the 
compression  and  reaction  phases  of  the  detonation  wave. 

The  ZND  model  is  an  oversimplification.  With  the  theory  of  liquids  suf¬ 
ficiently  advanced  that  quantitative  calculations,  including  nonideal  effects, 
can  be  made,  there  is  no  real  need  for  an  Inaccurate  treatment. 

Molecular  dynamics  simulations  have  shown  that  strong  dense-fluid  shock- 
waves  can  be  described  fairly  well  by  the  Navler-Stokes  equat1on(5).  Figure  2 


Indicates  the  relatively  good  agreement  between  the  Navler-Stokes  continuum 
calculations  and  the  atomistic  simulations.  The  principle  disagreements  are 
these:  The  wave  Is  about  30<  broader  than  the  continuum  predictions  and  the 
velocity  distribution  is  considerably  more  complex  than  is  the  equilibrium 
Maxwel l-Boltzmann  distribution.  The  alternative  Mott-Smith  approach(8),  in 
which  the  velocity  distribution  is  a  weighted  average  of  the  reactant  and 
product  Maxwellian  distributions,  grossly  overestimates  the  nonequilibrium 
character  of  the  wave. 

We  initially  planned  to  carry  out  a  similar  comparison  for  nitric  oxide,  a 
simple  liquid  which  undergoes  the  reaction 

2N0  =  N2  +  02.  (1) 

Exploratory  calculations  fitted  to  the  kinetics  of  the  reaction  showed  that  the 
reaction  zone  extends  for  thousands  of  angstroms,  far  too  long  for  molecular 
dynamics  simulation.  The  liquid  detonation  wave  profile  that  we  display  in 
Figure  3,  together  with  the  ZND  approximation  to  that  profile,  was  calculat¬ 
ed^)  with  an  activation  energy  approximately  half  that  of  nitric  oxide.  This 
substantially  reduced  the  spatial  extent  of  the  detonation  wave.  Under  this 
artificial  assumption  there  is  a  noticeable  quantitative  discrepancy  between 
the  simple  two-step  model  and  the  solution  of  the  coupled  equations.  There  is 
no  such  dramatic  difference  for  nitric  oxide,  at  least  for  a  one-dimensional 
wave. 

The  instability  of  one-dimensional  low-density  detonation  waves  is  well- 
established  experimentally(lO) .  Because  this  effect  is  at  least  as  important 
as  that  of  the  hydrodynamic  transport  coefficients  we  see  that  further  study  of 
the  one-dimensional  detonation  wave  for  nitric  oxide,  is  not  feasible,  using 
molecular  dynamics. 


Vibrational  energy  transfer  Is  very  slow  in  diatomic  molecules  because  the 
single  vibrational  mode  Is  well-separated  from  the  "lattice  modes."  The  situa¬ 
tion  is  less  clear  cut  for  polyatomic  molecules  with  many  degrees  of  freedom. 
Can  the  solid-phase  initiation  problem  be  understood  using  equilibrium  chemical 
kinetics  and  equilibrium  energy  surfaces?  Are  nonequilibrium  velocity  distri¬ 
butions  and  intramolecular  energy  transfer  significant?  The  answers  to  these 
questions  are  unknown  despite  considerable  interest  and  debate(ll).  We  are 

* 

concentrating  our  efforts  on  enhanced  understanding  of  solid-phase  Initiation, 
using  the  methods  outlined  in  the  following  sections. 


III.  COMPUTER  SIMULATIONS  AT  HIGH  RATES  OF  STRAIN 

To  explore  the  rapid  deformation  associated  with  shockwaves  in  condensed 
fluids  and  solids  we  have  developed  special  methods  for  solving  the  equations  of 
motion  of  atoms  in  a  rapidly-deforming  fluid  or  solid(12).  These  methods  were 
originally  developed  and  applied  to  the  simulation  of  viscous  flows(13)  and  were 
later  applied  to  solid-phase  plast1city{14) .  In  either  case  the  macroscopic 
flow  is  described  by  a  given  strain-rate  tensor  Vu,  where  u  is  the  hydrodynamic 
flow  velocity.  To  Incorporate  the  flow  velocity  into  the  equations  of  motion,  a 
perturbation,  the  tensor  product  of  Vu  and  Doll's  tensor  (the  dyadic  composed 
the  particles'  summed  coordinate-momentum  product  qp)  is  added  to  the 
Hamiltonian  function.  Equations  of  motion  result  which  include  contributions 
proportional  to  the  strain  rate: 
q  =  (p/m)  +  q- Vu  ; 

p*F-Vu-p;  (2) 

E  =  -  VP;VU  . 

Equations  (2)  describe  adiabatic  homogeneous  deformation.  This  scheme  has 
been  used  to  make  quantitative  estimates  of  the  shear  and  bulk  viscosities  for 
dense  fluids. 

Inhomogeneous  compressions  must  be  treated  differently.  The  system  is 
sheared  or  compressed  by  moving  periodic  images(5).  Either  method  could  be 
applied  to  chemically-reacting  solids.  We  expect  to  explore  the  homogeneous 
uniaxial  compression  of  two  systems:  The  simple  reaction  model  described  in  the 
next  section  and  a  "realistic"  three-dimensional  polyatomic  model  of  a  detonat¬ 
ing  solid.  This  latter  application  should  allow  us  to  map  out  the  path  followed 
by  energy  during  the  equilibration  process. 
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IV.  ATOMIC  DISTRIBUTION  FUNCTIONS  IN  SOLIDS 


We  are  currently  investigating  a  simple  model  for  reacting  atoms:  Two 
atoms  sufficiently  close  together  react,  releasing  an  energy  Q.  Equilibrium 
simulations  (with  Q*0)  have  been  carried  out  to  determine  the  density,  strain, 
and  temperature  dependence  of  the  "reaction  rate."  Figure  4  indicates  the 
substantial  dependence  on  strain  biaxiality.  Pastine's  model(15),  which 
contains  an  adjustable  collision  frequency  w,  can  be  used  to  fit  the  uniaxial 
compression  results. 

A  more  fundamental  approach  incorporates  the  distribution  of  neighboring 
atoms.  In  the  quasiharmonic  case  this  can  be  done  exactly.  We  have  developed 
a  numerical  method  for  the  partial  integration  of  the  Boltzmann  factor, 
exp(-<I|/kT) ,  over  all  but  two  particle  coordinates,  leading  to  an  atom-atom 
distribution  function.  This  function  is  a  product  of  Gaussians. 

Because  the  calculation  Is  somewhat  intricate  it  is  important  to  have  a 
check  available.  The  close-packed  lattice,  with  nearest-neighbor  interactions. 

Is  ideal  for  this  purpose.  Because  the  potential  energy  of  a  nearest-neighbor 
pair  In  a  D-dimensional  close-packed  crystal  is  kT/(D+l),  the  mean-squared 
deviation  in  the  nearest-neighbor  spacing  should  be  2kT/(D+lK,  where  <  is  the 
nearest-neighbor  force  constant.  We  have  verified  that  this  result  is  repro¬ 
duced  exactly  by  our  numerical  approach,  and  find  that  the  convergence  of 
higher-order  and  transverse  correlations  is  relatively  rapid  and  smooth  as  the 
number  of  particles  in  the  periodic  cell  is  increased.  In  the  two-dimensional 
case,  for  example,  the  ratios  of  the  nearest-neighbor  mean-squared  longitudinal  I 

to  transverse  fluctuations  for  9,  16,  25,  and  36-atom  crystals  are  0.667, 

0.688,  0.695,  and  0.696(16). 

We  expect  to  apply  this  same  technique  to  the  atom-atom  distribution  I 

function  for  polyatomic  crystals  undergoing  uniaxial  compression.  This  analog 
of  Pastine's  model* should  prove  useful  in  estimating  the  dependence  of  reaction 
rate  on  the  macroscopic  strain  rate.  _n-  j 

i 
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FIGURES: 


1.  Diatomic  reactant  and  product  Hugoniots  calculated  with  the  Lennard-Jones 
equation  of  state.  The  well  depth  and  collision  diameter  are  e  and  a, 
respectively.  The  Initial  volume  is  0.8442Na^.  The  initial  temperature 
is  0.722e/k. 

2.  Shockwave  profile  for  a  strong  fluid  shock  carrying  the  triple-point  argon 
to  a  final  temperature  of  about  one  electron  volt.  The  smooth  curves  are 
the  Navier-Stokes  solution.  The  points  correspond  to  data  from  nonequili¬ 
brium  molecular  dynamics. 

3.  Detonation  profile  calculated  using  (a)  the  hydrodynamic  equations.  Includ¬ 
ing  viscosity  and  heat  conduction,  and  (b)  using  the  ZND  model.  The  activa 
tation  energy  is  artifically  reduced  to  enhance  the  difference  between  the 
two  approaches.  The  heat  of  reaction  is  100c  and  the  activation  energy  is 
70e. 

X.  Reaction  rate  as  a  function  of  density  for  a  simple  binary-collision  model. 
Both  uniaxial  and  hydrostatic  compressions  are  shown  at  a  temperature  about 
half  the  melting  temperature.  Pastine's  model  for  the  rate  is  shown  as  the 
full  curve. 


Ij.  tU  K  iLtL'A'iI.'im  Vi 


F^uCt 


T)ifrToMlC  Ish^ARb-Ooi 
•frVToftfrTtoti  HuqoftiQT 

ti<rVv/0-  O.S^Z 
kT0/s  •-  o.?ll 


7 


-L\- 


